require(foreign)
require(MatchIt)
require(rbounds)
require(causalsens)
gender <- read.dta("genderjudging.dta")
gender$judge_vote <- 1 * (gender$judge_vote=="Liberal")
gender$gender_judge <- 1 * (gender$gender_judge == "Female")
gender$jud_experience <- 1 * (gender$jud_experience == "Experienced")
gender$minority_judge <- 1*(gender$minority_judge == "Minority")

gender$non_employ <- NULL
gender$hostile_work <- NULL

gender$treat <- 1 * (gender$num_fem > 0)

men <- subset(gender, gender_judge == 0)

## Propensity score model
pmod <- glm(treat ~ jcs + I(jcs^2) +minority_judge + minority_judge * jcs+
            as.factor(circuit) ,
            data = men, family = binomial())
men$pscores <- predict(pmod, type = "response")

## Matching model
m.out <- matchit(treat ~ minority_judge + jcs + jud_experience +
                   confirm_yr , data = men, distance = men$pscores,
                 exact = c("dec_year", "lower_dir"),
                 replace = TRUE, ratio = 6)
summary(m.out)
m.data <- match.data(m.out)

## Outcome model
y.out <- lm(judge_vote ~ treat, data = m.data, weights = weights)


## Rosenbaum SA
dis <- table(m.data$judge_vote, m.data$treat)
binarysens(x = 84, y=160, Gamma = 1.6, GammaInc = 0.05)

## Causal sensitivity analysis
alpha <- seq(-.5, 0.5, by = .01)
one.ate <- causalsens(y.out, pmod, ~ jcs + year_birth + minority_judge + lower_dir,
                       data = m.data, alpha = alpha)
one.att <- causalsens(y.out, pmod, ~ jcs + year_birth + minority_judge + lower_dir,
                       data = m.data, alpha = alpha, confound = one.sided.att)
align.ate <- causalsens(y.out, pmod, ~ jcs + year_birth + minority_judge + lower_dir,
                       data = m.data, alpha = alpha, confound = alignment)
align.att <- causalsens(y.out, pmod, ~ jcs + year_birth + minority_judge + lower_dir,
                       data = m.data, alpha = alpha, confound = alignment.att)

## Figure 2
par(mfrow = c(2,2))
plot(one.ate, bty = "n", main = "One-sided bias (ATE)", ylim = c(-0.4, 0.7))
plot(one.att, bty = "n", main = "One-sided bias (ATT)", ylim = c(-0.4, 0.7))
plot(align.ate, bty = "n", main = "Alignment bias (ATE)", ylim = c(-0.4, 0.7))
plot(align.att, bty = "n", main = "Alignment bias (ATT)", ylim = c(-0.4, 0.7))
